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ABSTRACT 


The  serru-Lagrangian,  semi-implicit  method  is  used  to  model  the  one  dimensional  shal¬ 
low  water  system  of  equations  with  surface  topography.  The  forecasts  are  compared  to 
finite  difference  and  semi-Lagrangian,  explicit  forecasts.  In  the  first  experiment,  a  non¬ 
rotating  system  is  considered.  The  semi-Lagrangian,  semi-implicit  model  agrees  very 
well  with  hydraulic  jump  theory,  while  the  semi-Lagrangian,  explicit  model  exhibits  ex¬ 
cessive  smoothing  and  the  finite  difference  model  breaks  down  when  the  nonlinear 
interactions  become  too  large.  In  the  second  experiment,  the  system  is  allowed  to  rotate 
to  examine  the  effect  of  rotation  on  the  formation  of  topographically  induced  hydraulic 
jumps.  Although  further  study  is  necessary,  it  is  clear  that  rotation  retards  the  devel¬ 
opment  of  the  low  pressure  to  the  lee  of  the  obstacle.  A  larger  domain  and  higher  spatial 
resolution  are  needed  for  more  detailed  simulation  of  hydraulic  jumps. 
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I.  INTRODUCTION 


There  are  many  numerical  methods  available  for  atmospheric  modelling,  but  two 
features  of  the  semi-Lagrangian,  semi-implicit  (SLSI)  method  seem  to  indicate  that  it  is 
particularly  well  suited  for  numerical  weather  prediction. 

The  first  of  these  traits  is  efficiency.  Because  the  semi-Lagrangian  method  approxi¬ 
mates  the  total  time  derivative  and  the  terms  responsible  for  the  high  frequency  gravity 
waves  are  treated  implicitly,  the  time  step  is  not  limited  b  the  Courant-Friedrichs-Lewy 
stability  criterion.  Thus,  a  relatively  large  time  step  can  be  used  with  the  SLSI  scheme 
to  produce  a  timely  forecast. 

Secondly,  a  recent  investigation  by  Kuo  and  Williams  (1989)  has  shown  that  while 
more  traditional  forecast  models  behave  poorly  when  large  gradients  develop  in  a  rela¬ 
tively  narrow  zone,  the  semi-Lagrangian  technique  can  still  produce  relatively  accurate 
solutions  when  a  scale-collapse  occurs.  This  is  an  important  attribute  for  many 
meteorological  applications,  such  as  predicting  the  formation  of  squall  lines  and  fronts, 
or  investigating  mountain  effects. 

In  this  study,  the  SLSI  method  is  used  to  model  the  shallow  water  system  of 
equations  with  topography  in  order  to  evaluate  the  handling  of  mountain  effects  by  the 
SLSI  method.  One  set  of  experiments  is  conducted  for  the  non-rotating  system.  Several 
sets  of  initial  conditions  are  used,  some  of  which  are  predicted  to  produce  hydraulic 
jumps.  The  SLSI  forecasts  are  compared  with  the  forecasts  of  two  other  models:  a 
leapfrog  finite  difference  scheme  (FDEX),  and  a  semi-Lagrangian.  explicit  (SLEX) 
model.  Next,  the  experiments  are  repeated  for  a  rotating  system  to  examine  the  elfects 
of  rotation  on  the  formation  of  orographically  produced  jumps.  It  is  known  that  rota¬ 
tion  can  prevent  the  development  of  hydraulic  jumps  in  a  fiat  bottom  system  (Williams 
and  Hori,  1970). 

In  the  next  chapter,  the  development  of  the  semi-Lagrangian  technique  is  reviewed, 
along  with  theory  on  the  formation  of  hydraulic  jumps.  In  Chapter  III,  the  governing 
equations  and  formulation  of  the  FDEX,  SLEX,  and  SLSI  models  are  discussed.  The 
initial  conditions  and  results  of  the  experiments  are  presented  in  Chapter  IV,  and  con¬ 
clusions  and  recommendation  for  further  study  are  in  Chapter  V. 
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II.  BACKGROUND 


A.  THE  SEMI-LAGRANGIAN  SCHEME 

In  a  pure  Lagrangian  scheme,  such  as  Fjortoft  (1952)  proposed,  one  set  of  fluid  el¬ 
ements  is  tracked  for  the  entire  integration  period.  Thus,  a  parameter,  Q,  that  is  con¬ 
served  following  the  fluid  such  that 


80  SO 

“■+  ys-r  =0, 


ot 


os 


(2-1) 


where  s  is  the  direction  of  the  flow  and  V  is  the  velocity,  remains  constant  with  inte¬ 
gration.  With  the  more  common  Eulerian  methods,  a  different  set  of  fluid  elements  is 
evaluated  at  each  time  step,  representing  a  different  distribution  of  the  advective  pa¬ 
rameters.  Fjortoft's  goal  was  to  develop  a  scheme  in  which  a  large  time  increment,  Ar, 
could  be  used.  In  regions  of  strong  wind  shear  or  after  long  integration  periods,  however, 
the  initial  grid  can  become  greatly  distorted,  as  illustrated  in  Fig.  1.  Data  points  may 
cluster  in  a  relatively  small  area,  leaving  large  areas  unanalyzed. 


To  avoid  the  problem  of  distortion  and  still  make  explicit  use  of  equation  (2-1), 
Wiin-Nielsen  (1959)  proposed  a  quasi-Lagrangian,  or  trajectory,  scheme.  The  fluid's 
properties  are  evaluated  on  a  regular,  stationary  grid  at  each  time  step  by  determining 
the  off-grid  positions  from  which  the  fluid  elements  originated. 

The  semi-Lagrangian  method,  first  proposed  by  Sawyer  (1963),  is  a  modified  version 
of  the  Wiin-Nielsen  scheme.  In  this  study,  a  one-dimensional,  three  time  level  version 
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of  the  semi-Lagrangian  method,  similar  to  the  two-dimensional  scheme  described  by 
Robert  ( 19S1),  is  used.  For  each  parameter,  Q,  such  that 


dQ  cO 

7r+(/ir=* 


(2-2) 


where  R  represents  all  the  remaining  terms,  the  left  hand  side  is  approximated  at  the  grid 
point,  .r,  ,  by 


Q{xht  +  At) 


Qjxi  -2 aht  -  At) 
2  At 


=  R, 


(2-3) 


where,  by  the  midpoint  rule, 


ai  =  *t  •  Ufa  —ctrt).  (2  -  4) 

The  o' s  must  be  solved  for  iteratively.  For  each  n, 

a*  =  At  •  U(x;  —  a?~\t).  (2  —  5) 


Subtracting  equation  (2-4)  from  (2-5)  and  applying  the  Lipschitz  condition  leads  to 


Ik-"1  -kl 


=  At  ■ 


cU 

cx 


(2-6) 


Thus,  equation  (2-5)  will  converge  whenever  the  right  hand  side  of  Equation  (2-6)  is  less 
than  one.  Kuo  and  Williams  (19S9)  argue  that  no  more  than  three  iterations  are  needed, 
if  this  condition  is  met,  since  equation  (2-4)  is  only  third  order  accurate  with  respect  to 
At  and  equation  (2-5)  becomes  order  At  more  accurate  with  each  iteration.  Pudykiewicz 
and  Staniforth  (19S4)  point  out  that  the  maximum  horizontal  wind  shear  in  most 
meteorological  models  is  of  the  order  lO^-1,  so  a  time  step  of  At<2  hours  is  necessary 
for  convergence.  By  contrast,  the  time  step  in  most  Eulerian  models  used  in 
meteorological  applications  is  restricted  by  the  Courant-Friedrichs-Lewy  (CFL)  stability 
criterion, 


C  = 


u 


At 

Ax 


<  1, 


(2-7) 
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where  C  is  the  Courant  number.  The  time  step  is  limited  by  both  the  maximum  wind 
speed,  Cmax,  and  the  resolution  of  the  model.  For  a  given  l'max,  the  higher  the  resolution, 
the  smaller  the  time  step. 

B.  DEVELOPMENT  OF  HYDRAULIC  JUMPS 

In  this  study,  the  one-dimensional  shallow  water  system  of  equations  is  used  to  ex¬ 
amine  the  nature  of  the  flow  of  an  incompressible,  hydrostatic,  homogeneous,  inviscid 
fluid  over  an  obstacle  in  both  a  non-rotating  and  rotating  system.  Under  certain  con¬ 
ditions,  it  is  possible  hydraulic  jumps  may  form  in  one  or  both  cases. 

1.  The  Non- Rotating  S  stem 

Houghton  and  Kasahara  (1968)  studied  the  non-rotating  shallow  water  system, 
illustrated  in  Fig.  2,  and  determined  the  conditions  which  would  lead  to  the  formation 
of  hydraulic  jumps.  When  there  is  no  Coriolis  force,  i.e./  =  0,  then  the  one-dimensional 
shallow  water  equations  are 


cu  cu  ch  A 

-T-+U  —  +  g^^  =  0 

Cl  cx  cx 


(2-8) 


and 


ell 

ci 


ell  ,,  cu  A 
-t-  U— —  +  H  ~  =  0, 

cx  cx 


where  u  is  the  horizontal  velocity,  H  is  the  depth  of  the  fluid  and 

h'  =  H  —  hB. 


(2  -  9) 


(2-10) 


t—  represents  the  height  perturbation  on  the  free  surface. 

C  X 

At  t  <  0,  the  fluid  is  at  rest,  and  the  free  surface  is  flat.  At  i  —  0,  the  fluid  is 
given  a  constant  velocity  of  u0  throughout  the  entire  domain.  At  t  =  oo,  the  fluid  is  at 
its  new  steady  state  given  by 

2  2 

^f-  +  g(Hss  +  hB)  =  ^-  +  gh  =  C]  (2-11) 

and 

ussHss  =  u0h  =  C2.  (2-12) 
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Eliminating  //„  from  equation  (2-11)  using  Equation  (2-12)  and  introducing  the 
following  dimensionless  parameters, 

1 1  1 1  p  it 

Fr  =  ==- ,  *  =  and  U.  =  -jj-  ,  (2-13) 

vV*  h 

into  the  resulting  equation  leads  to 

r  2  r  2 

-~-U^  +  U.{R 1)+  1  =0.  (2-14) 

For  a  given  Froude  number,  Fr,  R  can  be  plotted  as  a  function  of  U.  using  equation 
(2-14).  A  plot  of  F  versus  U.,  shown  in  Fig.  3  for  Fr2  =  0.1,  illustrates  there  are  three 
real  roots  to  equation  (2-14)  if  R  <  RCritical-  Only  one  °f  these  roots  is  physically 
meaningful.  If  R  >  RCriticau  only  one  solution  to  equation  (2-14)  exists,  but  since  this 
solution  has  no  physical  meaning,  a  hydraulic  jump  is  anticipated.  RCritical  can  be  ex¬ 
pressed  as  a  function  of  the  Froude  number: 

^■CRITICAL  =  ~2  ~  ~2  3  +  1’  (2—15) 
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A  plot  of  RCritical  versus  Fr,  illustrated  in  Fig.  4  indicates  the  three  classes  of  solution 
to  equation  (2-14)  identified  by  Houghton  and  Kasahara. 

Domain  I  is  the  suba  ideal  range  in  which  R  <  Rcritical  and  Fr  <  1.  The  steady- 
state  free  surface  height  of  a  fluid  flow  which  meets  these  criteria  will  dip  over  the  ob¬ 
stacle.  The  velocity  will  increase  over  the  obstacle  but  will  remain  less  than  u.,  the  speed 
that  corresponds  to  the  condition  R  =  RCRincAi •  Domain  II  is  the  'jump  region'  in  which 
R  >  Rcritical  and  a  hydraulic  jump  forms.  Domain  111  is  the  supercritical  range  in  which 
R  <  Rcritical >  but  Fr  >  1.  At  its  steady  state,  the  free  surface  height  rises  over  the  ob¬ 
stacle,  and  the  velocity  decreases  but  remains  greater  than  u.. 
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Figure  4.  R critical  versus  Fr 

2.  The  Effect  of  Rotation 

If  the  coordinate  system  is  rotating,  i.e.  / ^  0,  then  the  shallow  water  system 
equations  for  a  flat  lower  boundary  and  no>»  dependence  are 

i7+“ir+*  i7->=0’  <2-'6> 

f  +  U£+/u-0  (2-. 7) 

and 
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(2-18) 


cH 

8t 


^  dH 
-b  u  ~  + 

C  X 


H-rr- 

cx 


=  0. 


If  the  variables  are  scaled  as  follows: 


x  =  Lx,,  t  = 


yjghm 


u  =  Uu, 


(  Fr  \ 


H  =  h  +  U  Yhs - 


(2-19) 


n  u  J  r  u 

R„  =  ~  and  F  =  — ==- , 


0  fL 


V 


gh 


then  equations  (2-16).  (2-17)  and  (2-1S)  become 


cu,  8u,  ch, 

■+  Fus-rL  +  r 


01 ' 


ox. 


*  ex. 


■4-o. 


(2  -  20) 


cv 


o  I'p  C  lr 

— —  +  Fus  — —  +  =  0 


cr 


ox. 


(2-21) 


and 


c/t,  c/tr  cu,  cu, 

s  +F(us-^-  +  hs-rrL)  +  -^-  =  0. 


crr 


cjrc 


c.t,  c.v. 


(2  -  22) 


Unless  the  Coriolis  term  dominates,  hydraulic  jumps  may  form  due  to  the  nonlinear 
advection  terms.  The  no-jump  condition  can  be  expressed  as 


F 

rI 


(2  -  23) 


or 


F  =  A  •  Rq, 


(2-24) 


where  A  is  some  number  much  larger  than  one.  The  Williams  and  Hori  (1970)  investi¬ 
gation  found  that  A  is  in  the  range  of  6.5  to  7.0. 
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3.  Scale  Collapse  Problem 

Computational  difficulties  may  arise  when  a  numerical  model  must  represent  a 
physical  process  in  which  a  large  gradient  or  discontinuity  forms.  Kuo  and  Williams 
(19S9)  considered  a  simple  scale-collapse  case,  for  which  the  governing  equation  was  the 
one  dimensional  nonlinear  advection  equation 


cu 

ct 


,  A 

+  u  ~  =  0, 

cx 


(2-25) 


and  the  initial  condition 

«W,=o  =  ~  tan-1  {x  -  x0).  (2  -  26) 

This  problem  can  be  solved  for  analytically,  and  that  solution  leads  to 

cu  ,  1  , 

~t —  r_r  = - : - ►  —  oo  as  t  —+  l  .  U  —  2/ ) 

cx  x~x»  /  —  1  v 

Kuo  and  Williams  compared  the  second-order  centered  finite  difference,  the 
Chebyshev-Tau  and  the  semi-Lagrangian  solutions  of  equation  (2-25)  and  (2-26).  The 
scale-collapse  at  ,r  =  produces  large,  fluctuating  errors  throughout  the  entire  domain 
when  either  the  finite  difference  or  Chebyshev-Tau  method  is  used.  The  semi- 
Lagrangian  method  responds  differently.  As  t  approaches  one  from  below, 

AMI  4^  ||  >1,  (2  —  28) 

cx 

for  Therefore,  equation  (2-5)  will  not  converge  in  this  region.  However,  if  a 

fixed  number  of  iterations  {n  =  3)  is  used  to  solve  for  the  a,' s,  the  errors  which  result 
when  and  where  equation  (2-28)  is  true  remain  near  the  scale  collapse  zone.  Thus,  the 
semi-Lagrangian  method  is  still  useful  even  after  discontinuities  develop. 

In  this  study,  two  semi-Lagrangian  schemes  for  the  shallow  water  equations  are 
compared  with  the  traditional  leap-frog  finite  difference  method.  In  the  first,  the  R  terms 
in  equation  (2-2)  are  treated  explicitly;  that  is,  they  are  evaluated  at  (x,—  a„t).  This 
scheme  is  expected  to  handle  large  shear  zones,  but  the  freedom  from  the  CFL  criterion 
cannot  be  exploited,  because  the  time  step  is  still  limited  by  the  gravity  wave  terms.  The 
second  scheme  is  a  semi-Lagrangian,  semi-implicit  scheme  similar  to  the  one  presented 
by  Pudykiewicz  and  Staniforth  (19S4),  but  modified  to  include  topographical  effects. 
The  terms  responsible  for  the  gravity  waves  are  treated  implicitly  by  averaging  between 
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(x„t  +  A i)  and  (jcf  -2 a„t  -  A/).  The  remaining  terms  are  evaluated  at  (x,  -a„i).  The 
semi-Lagrangian,  semi-implicit  scheme  requires  more  computational  effort,  but  it  per¬ 
mits  a  much  greater  time  step. 
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III.  FORMULATION  OF  THE  MODELS 


A.  THE  GOVERNING  EQUATIONS 

The  governing  equations  for  the  shallow  water  model,  illustrated  in  Fig.  2,  are 


cu  ,  cu  , 

-7-  +  u—  +  g 

ct  cx 


cH_  M B 
ox  dx 


—fv  =  0, 


(3-1) 


CV  ,  CV  ,  / 

-r-  4-  u  — - 1-  g 

ct  cx  ' 


IK  +  S^)+/U  =  o 

cy  cy  I 


(3-2) 


and 


^¥-  +  -?-(uH)  +  4-{vH)  =  0. 

ct  cx  cy 


(3-3) 


In  the  one-dimensional  model,  the  mean  flow,  1'0,  is  zonal  and  in  geostrophic  balance 
with  the  mean  free  surface  height,  such  that 


u 


O 


.e  rh 
f  cy 


(3-4) 


and 


v 


O 


S  Ch 

f  cx 


=  0, 


(3-5) 


and  perturbations  in  the  dependent  variables,  u',  v'  and  h'  are  in  the  ar  direction  only. 
Multiplying  equation  (3-3)  by  g  and  substituting 

u(x,t)  =  u0  +  u’(x,t),  (3  -  6) 

v(x,t)  =  /(ar,/)  (3-7) 


and 


4>(x,t)  =  4>{y)  -  <t>B(xy)  +  4>'(x,r)  =  gH  (3  -  8) 

into  equations  (3-1),  (3-2)  and  (3-3)  yields  the  following  set  of  equations 
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(3-9) 


CU  ,  /  M  ru  C(j) 

ci  K  0  c.v  c.t  ' 


^  +  (uo  +  »')4L-+/«'  =  0 
Cl  cx 


and 


c<t>' 


ci 


+  (u0  +  u') 


cd)' 


cx 


c4>b 


CX 


+  (<f>  —  4>  B  +  4> 


i%) 


+ 


(3  -  10) 


(3-11) 


Tor  convenience,  the  bottom  topography  is  assumed  to  be  parallel  to  the  mean  free 
surface  height  in  they  direction,  so  that 


rh  ?hB 


(3-12) 


and  equation  (3-11)  reduces  to 


54)’  (  c&’ 

■  +  K  +  «')( 


ct 


cx 


— f—  )  +  (<■>-  (pB  -f  <£')  -TT-  =  0.  (3-13) 


cx 


cx 


B.  FINITE  DIFFERENCE  SCHEME 

One  of  the  simplest  numerical  methods  is  finite  differencing  centered  in  time  and 
space,  the  leapfrog  scheme.  On  the  staggered  grid,  equations  (3-9),  (3-10)  and  (3-13) 
approximated  by  this  method  become 


u'(x.i  -I-  Ai)  -  ''(x.i  —  At) 

_ 


+  (u0  +  u'(x,t)) 


u'(x  +  Ax.t)  —  u'(x  —  Ax.t) 
2Ax 


Ax 


+f 


v'(x  -  ~  j)  +  v'(x  +  -y ,/) 


.(3-14) 


v'fv.r  +  Ar)  —  v'(x,r  —  Ar) 
2At 


+  u(x,t) 


v'(x  +  Ax.t)  —  v'(x  —  A.v.r) 
2Ax 
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r  u'(*  + 


(3-15) 


Ax  s  ,,  Ax  , 
—  J)  +  u{. x —  ,r) 


and 


4>'{x,t  +  A/)  -  ct>'{x.i  -  A/) 
- - +  u{x,t) 


4)'(x  +  Ax,r)  —  (j)'(x  —  Ax.r) 
2Ax 


where 


u(x,t) 


<t>B{x  +  Ax)  -  4>b{x  -  Ax) 
2  Ax 


-  (<5  -  <bB{x)  +  <t>'(x,t)) 


,,  .  Ax  ^  ,,  Ax  s 

u  (x  +  —  ,r)  -  u  (x  -  —  ,i) 

Ax 


(3-  16) 


u(x.t)  =  u0 


u'(x  -  ,r)  +  u'(x  +  A—  j) 


(3-17) 


C.  SEMI-LAGRANGIAN,  EXPLICIT  SCHEME 

The  second  numerical  scheme  uses  the  semi-Lagrangian  method  described  by  Robert 
(1981),  modified  for  one  dimensional  How  over  a  surface  with  topography.  The  re¬ 
maining  terms  are  centered  in  space  about  (x  —  a,i).  Equations  (3-9),  (3-10)  and  (3-13) 
approximated  by  this  method  become 


u'(x.i  +  Ar)  —  u'(x  —  2a. t  —  A  r) 
2Ar 


S(j)' 

ox 


l(x-a.r)  +>'(* 


a,t), 


(3-  U>) 


v'(x,r  +  Ar)  —  v'(x  —  2  a,t  —  Ar) 

2A 7 


=  -fu'(x  -  a,i ) 


(3-19) 


and 

<t>'{x,t  +  Ar)  —  <t>'(x  —  2a, r  —  Ar)  4>b(x)  ~  <£b(x  —  2a) 
2At  2Ar 
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(3  -  20) 


where 


6Q_  SO  Qi*  +  —  )-Q(x - 2~) 

cx  <Lv  ^  Ax 


(3-21) 


C(y  Q\X 

Equation  (3-21)  is  used  to  estimate  and  at  the  grid  points,  and  then  an  in¬ 
terpolation  scheme  is  used  to  evaluate  these  terms  at  (jc  —  a).  An  alternative  m  .hod 
would  be  to  approximate  the  derivatives  immediately  at  (jr  —  a)  with  the  following 
equation: 


io  fi(*  -  a  +  )  -  Q(x  -  a  -  ) 

77  '<*-«> - ^77 - —  •  <3-2:» 

Although  this  method  appears  more  direct  than  using  equation  (3-21),  it  requires  twice 
as  many  interpolations  and  therefore  more  computational  effort. 

Note, 


_  c4>b  c4>b  _  c<t>B  4>b(x)  -  $B(x  -  2d) 

u  — —  =  — —  -I-  u  — —  2: - — - 

c.r  ct  cx  2  At 


(3-23) 


since  there  is  no  local  time 


change  in  the  bottom  topography,  i.e., 


c4>b  . 

— —  is  zero. 
ct 


D.  SEMI-LAGRANGIAN,  SEMI-IMPLICIT  SCHEME 

This  scheme  is  similar  to  the  semi-Lagrangian,  explicit  scheme  presented  above  ex¬ 
cept  the  gravity  wave  terms  are  treated  semi-implicitly.  The  development  of  this  scheme 
is  similar  to  the  development  of  the  semi-Lagrangian,  semi-implicit  scheme  developed 
by  Staniforth  and  Temperton  (1986)  for  a  flat  bottom  system.  The  finite  difference 
equations  now  become 


u'(x,r  +  At)  —  u'(x  —  2a,  t  —  At) 
2At 


1 

2 


S4>'  Scf)' 

Sx  w+ao+  sx 


+fv'(x  -  a,t), 


v'(jr,r  -f  At)  —  v'(jc  —  2a, t  —  At) 
2Ax 


-fu'{x  -  a,t) 


(3  -  24) 


(3-25) 
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and 


d>'(x.t  +  At)  —  d>'(x  —  2a,t  -  At)  4) B(x)  -  d> B(x  -  2a) 

2 At  2At  =  “  * 

-  y  [(*  -  ***))  77 !(,.,«,)  +  (*  -  -  *0)  x-  <3  -  26> 

The  divergence  term  at  (t  +  At)  is  eliminated  from  equation  (3-26)  by  first  finite  differ¬ 
entiating  equation  (3-24)  with  respect  to  x,  which  yields, 

Su’  |  _  dtt'  .  _  A 

$x  W+Ar)  —  t^r  l(jc-2fl,r-Af)  &1 


o  4>’  , 

.  2  W+Ar) 

dx 


+ 


dx 2 


At) 


+  dx 


(3  -  27) 


where 


d2(t)'  4)'{x  4-  A.r)  —  2</>'(x)  4-  4>'(x  —  A.r) 

<5.r2  A.t2 


(3  -  28) 


then  substituting  equation  (3-27)  into  (3-26)  and  collecting  all  (.r.t-l-At)  terms  on  the 
left-hand  side 


ci2  _  1 

6x'  A r{4>  -  4>b(x)) 


)0'(.r.t  4-  At)  =  — 


•  + 


dx 


Arid)  -  d>B[x)) 


((^'(a-  —  2a,t  —  At) 


4>b(x)  ~  4>b(x  -  2«)  /  1  ^  “  2a)  A  Su'  , 

A  t\t-4)B(x))  +\^+  JSt^'-^ 


2<l>'(x-a)  Si/ 
A t(0  -  <f>B(x))  Sx 


I  (x-a,t) 


+  2f 


Sv'  | 

(5  V  kx-a.l)' 


(3  -  29) 


or,  expanding, 


A.r 


4>'(x  +  Ax)  -  2  +  — r== 

V  At2(<£  -  4>b{. r)) 


</>'(*)  +  0'(*  -  Ax) 


Ax‘ 


If+Af 
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0'(.r  —  2a  +  Ax)  —12  — 


Ax 


Af(<t>  -4>b(x-  2a)) 


]C>'(.y  —  2d)  +  4>'(x  —  2a  —  A.y) 


Ax 


-  <M-r~ 2g)  |  f  i  <t>-  -  ?a) 

Ai7{(j>  -  y))  V  At  A[($  -  4>b(x)) 


u'(x  —  2a  + 


Ax  x  _ 


)  -  u'(x  -2a-  -r- ) 


Ax 


2 (£'(•*•-  a) 

A t{4>  -  <M*)) 


«'(.*  —  a  + 


A.y 

2 


)  —  u'(x  —  a  — 


Ax 


+f 


v'(.y  +  a  +  A.y)  —  v'(.y  —  a  —  A.y) 
A.y 


(3  -  30) 


Equation  (3-30)  can  be  written  in  matrix  form  as 

| >*]{</>',}  =  {£,},  (3-31) 

where  { B ,}  are  all  the  term  ;  on  the  right-hand  side  of  equation  (3-28)  and  the  matrix  A 
is 


where 


a,  1  0...  ...0  1 

1  a2  1  0...  ...0 


M  = 


0...  ...0  1  a  v_,  1 

1  0...  ...0  1  a  v 


_  Ay2  \ 

A r(<t>  -  A*))  ) 


To  solve  equation  (3-31),  it  is  necessary’  to  invert  the  matrix  A  : 


{<*>';}  =  [XT1 


(3  -32) 


(3-33) 


(3  -  34) 
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However,  since  matrix  A  has  no  time  dependent  elements,  the  inverse,  [/[]-',  must  be 
solved  for  only  once. 

E.  THE  SPATIAL  GRID 

u' ,  v'  and  <£',  are  evaluated  on  the  staggered  grid,  scheme  C  of  Arakawa  and  Lamb 
(1977),  illustrated  in  Fig.  5.  The  use  of  scheme  C  further  limits  the  maximum  time  step 
for  the  explicit  models,  already  restricted  by  the  CFL  criterion,  by  a  factor  of  two,  as 
Haltiner  and  Williams  (1980)  demonstrated  with  the  leapfrog  finite  difference  scheme. 
However,  they  also  showed  that  scheme  C  simulates  he  shorter  waves  far  better  than 
the  unstaggered  grid,  an  important  feature  when  investigating  orographic  effects.  Cyclic 
boundary  conditions  are  used  to  avoid  the  difficulties  associated  with  the  advection  of 
fluid  elements  from  outside  the  domain  onto  the  grid. 

An  interpolation  scheme  is  needed  to  evaluate  the  dependent  variables  at  {x  —  a)  and 
(x  —  2a).  but  the  resultant  smoothing  of  the  variables  between  grid  points  can  represent 
significant  errors.  Kuo  and  Williams  (1989)  have  shown  that  in  regions  where  the 
Courant  number  is  less  than  one.  interpolation  errors  dominate  time  truncation  errors. 
In  these  situations,  the  forecasts  become  more  accurate  as  A t  is  increased,  because  fewer 
time  steps  require  fewer  interpolations. 

In  this  study,  a  cubic  spline  interpolator  is  used.  Although  computationally  expen¬ 
sive,  the  cubic  spline,  because  of  it  accuracy,  has  been  a  popular  choice  for  use  with  the 
semi-Lagrangian  method  (Robert  (19S1),  Pudykiewicz  and  Staniforth  (19S4),  Ritchie 
(1987).  Kuo  and  Williams  (19S9).)  The  global  nature  of  the  cubic  spline  (  all  grid  values 
are  used  to  estimate  to  a  single  point),  led  Pudykiewicz  and  Staniforth  (19S4)  to  suggest 
that  this  scheme  would  be  inappropriate  where  locally  steep  gradients  occur  and  to 
speculate  that  a  local  interpolating  function  should  be  used  in  these  situations  to  mini¬ 
mize  the  Gibbs'  phenomenon.  However,  Kuo's  and  Williams'  (1989)  investigation  of  a 
scale  collapse  model  found  that  the  deviations  in  semi-Lagrangian  forecast  from  the 
analytical  solution  remained  small  and  confined  to  the  region  of  the  large  wind  conver¬ 
gence,  even  though  a  cubic  spline  interpolator  was  used. 
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Figure  5.  Staggered  Grid 
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IV.  EXPERIMENTS  AND  RESULTS 


In  the  first  experiment,  the  three  models'  forecasts  for  How  over  a  topographical 
ridge  in  a  non-rotating  system  are  compared.  The  same  cases  that  Petroliagis  (1988) 
investigated  with  a  Galerkin  model  are  examined  here.  All  three  domains  identified  by 
Houghton  and  Kasahara  (1968)  are  represented. 

In  the  second  experiment,  the  effect  of  rotation  upon  topographically  induced  jump 
formations  is  examined.  The  three  todels  are  run  with  the  same  sets  of  parameters  used 
in  the  first  experiment,  except  /=  0.000103  sec-1,  corresponding  to  the  Coriolis  force  at 
-15°  N. 

A.  INITIAL  CONDITIONS 

The  domain  of  integration,  L,  is  7276  km  in  the  x  direction,  divided  into  48  incre¬ 
ments  of  length  u  is  evaluated  on  the  odd  numbered  grid  points;  v  and  </>  are  eval¬ 
uated  on  the  even  numbered  grid  points.  The  boundary  conditions  are  cyclic,  such  that 

i/(0)  =  u(L),  v(0)  =  v(L)  and  0(0)  =  0(L).  (5-1) 


The  ridge  profile,  illustrated  in  Fig.  6  is  described  by 


"1 


4/ -20)  \ 
12  ) 


0 


for  20A.v  <  x  <  32Ajc 

elsewhere. 


(5-2) 
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The  five  cases  examined,  summarized  in  Table  1,  are  parameterized  by  the  mean 
height,  h  ,  the  ratio  of  the  mountain  peak  to  the  mean  height,  R,  and  the  Froude  num¬ 
ber,  Fr.  These  are  the  same  sets  of  parameters  that  Petroliagis  (1988)  used.  However, 
for  these  experiments  L  and  Ax  are  increased  tenfold  which  increases  the  time  scale,  so 
that  f  could  be  important  in  the  second  set  of  experiment. 


Table  1.  PARAMETERS  FOR  CASES  I  THROUGH  V 


CASE 

h  (m) 

R 

Fr 

u7  (m  s) 

D 

I 

1000. 

■0C 

0/ 

19.8 

I 

II 

1000. 

0.2 

0.4 

39.6 

I 

III 

1000. 

0.7 

0.3 

29.7 

II 

IV 

500. 

0.2 

0.8 

56.0 

II 

V 

400. 

0.05 

1.4 

87.7 

III 

9  t  ^ 

A  five-minute  time  step,  AtEX,  is  used  with  the  Finite  Difference  Explicit  (FDEX)  and 
Semi-Lagrangian  Explicit  (SLEX)  models.  A  one-hour  time  step,  AtSI,  is  used  with  the 
Semi-Lagrangian  Semi-Implicit  (SLSI)  model.  The  Semi-Lagrangian  Semi-Implicit 
model  is  also  run  at  a  higher  resolution  (HRSL)  for  comparison  with  the  other  models. 
For  HRSL,  Ax  is  reduced  by  a  factor  of  four,  and  AiSI  is  used. 

AtEX  is  limited  by  the  CFL  stability  criterion  based  on  the  speed  of  the  explicitly 
treated  external  gravity  waves.  The  mean  speed  of  these  waves,  c,  is  given  by 

c  =  |  u0  +  yjj)  | , 


=  l(Fr±l)V>  I- 


(4-3) 


A  mean  Courant  number  for  the  explicit  cases,  can  be  defined  as 


Cex  = 


2 


cJ±1ex 

Ax 


<  1. 


(4-  j 


As  discussed  in  chapter  3,  the  factor  of  two  arises  because  scheme  C  is  used.  CEX  is  ap¬ 
proximately  0.25-0.30  for  the  cases  considered. 

B.  RESULTS  I 

The  u  field  for  case  I  after  two,  four,  six  and  eight  hours  of  integration  are  shown 
in  Fig.  7.  All  the  models  produce  a  stationary  speed  maximum  centered  over  the  ridge, 
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which  the  theory  predicts  for  the  subcritical  range  (Domain  I).  FDEX  and  SLEX 
slightly  underdevelop  the  speed  maximum  while  SLSI  is  in  good  agreement  with  HRSL. 
All  the  models  show  the  development  of  a  secondary  maximum  that  propagates  east¬ 
ward.  After  eight  hours  the  downstream  disturbances  have  propagated  around  to  the 
east  side  of  the  ridge  and  begin  to  interfere  with  the  upstream  fields. 

As  anticipated,  a  low  pressure  region  develops  on  the  lee  side  of  the  mountain 
(Fig.  8).  The  high  pressure  field  to  the  east  of  the  ridge  remains  fairly  stationary  in  the 
SLSI  and  HRSL  forecasts.  This  feature  drifts  upstream  perceptibly  in  the  FDEX  and 
SLEX  forecasts  Theory  predicts  the  formation  of  a  wave  train  on  the  lee  side  of  the 
mountain.  FDEX,  SLEX  and  SLSI  all  lag  behind  the  high  resolution  model  in  devel¬ 
oping  these  waves.  SLSI  is  only  slightly  behind.  The  FDEX  becomes  appreciably 
damped,  but  SLEX  has  the  poorest  results  due  to  excessive  smoothing  after  numerous 
interpolations. 

Case  II  is  also  in  Domain  I,  however  the  mean  flow,  UM£ASt  is  stronger  and  the 
mountain  peak,  Hs.  is  higher.  Although  the  perturbations  are  relatively  larger,  the  u  and 
4>'  fields,  shown  in  Fig.  9  and  Fig.  10  for  two,  four,  six  and  eight  hours,  are  quite  similar 
to  their  counterparts  in  case  I. 

Case  III  and  case  IV  lie  in  Domain  II,  so  hydraulic  jumps  are  anticipated.  Fig.  11 
and  Fig,  12  illustrate  the  u  and  4>'  fields  at  three,  six,  nine  and  twelve  hours  for  case  III. 
All  models  forecast  a  wind  speed  maximum  over  the  ridge.  The  FDEX  solution  does 
not  appear  in  Fig.  1 1 C  and  Fig.  1 1 D  and  in  Fig.  12C  and  Fig.  12D,  because  FDEX  is 
unstable,  failing  to  converge  after  six  hours,  due  to  the  large  nonlinear  interactions  in  the 
vicinity  where  the  jump  is  developing.  The  HRSL  model  continues  to  deepen  and 
broaden  the  pressure  minimum  on  the  lee  side  of  the  ridge  until,  after  eleven  hours  of 
integration,  the  pressure  has  dropped  to  zero.  Once  the  free  surface  height  hits  the  bot¬ 
tom  topography,  the  integration  is  stopped.  The  perturbation  fields  develop  slightly 
slower  in  the  SLSI  model;  the  free  surface  hits  the  bottom  after  13  hours  of  integration. 
Initially,  the  SLEX  forecasts  are  very  similar  to  SLSI's.  With  time,  however,  the  SLEX 
perturbations  dampen  out,  again  due  to  interpolation  error. 

The  u  and  </>'  fields  at  three,  six,  nine  and  twelve  hours  for  case  IV  are  shown  in 
Fig.  13  and  Fig.  14.  The  amplitude  of  disturbances  are  increasing,  but  because  of  the 
cyclic  boundary  conditions,  there  is  interference  from  the  transient  part  of  the  solution 
before  a  jump  can  form. 

Case  V  is  in  Domain  III,  the  supercritical  range.  The  u  and  <p'  fields  from  all  the 
models  are  compared  in  Fig.  15  and  Fig.  16.  SLSI  again  agrees  well  with  HRSL.  Both 
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produce  the  speed  minimum  over  the  ridge  as  the  theory  predicts.  FDEX  and  SLEX 
develop  this  minimum  slightly  upstream.  All  the  models  show  a  pressure  minimum  to 
the  lee  of  the  obstacle.  The  SLSI  and  HRSL  remain  rather  stable,  until  the  transient 
solutions  propagating  to  the  east  circle  back  around  to  the  west  of  the  ridge.  The  SLEX 
model  exhibits  exceesive  smoothing  from  accumulating  interpolation  errors,  as  it  does 
for  all  the  other  cases.  The  FDEX  solution  becomes  unstable  in  this  case,  unable  to 
handle  the  nonlinear  interactions  when  Fr  exceeds  one. 

C.  RESULTS  II 

W  en  the  system  is  allowed  to  rotate,  perturbations  arise  in  v  as  well  as  in  u  and  4>. 
These  disturbances  in  v  are  out  of  phase  with  the  u'  field,  as  illustrated  in  Fig.  17  for 
case  I.  After  several  hours  of  integration,  perturbations  in  u  are  noticably  smaller  and 
the  pressure  drop  on  the  lee  side  of  the  obstacle  appears  not  as  deep  when  compared  to 
the  non-rotating  cases,  as  in  Fig.  18A  and  Fig.  18B  for  the  HRSL  results  for  case  II. 
Although  the  Coriolis  force  is  deflecting  some  energy  into  they’  direction,  cases  III  and 
IV  are  still  unstable.  For  both  cases,  the  amplitudes  of  the  disturbances  in  the  vicinity 
of  the  ridge  continually  increase.  This  is  illustrated  in  Fig.  19  for  the  total  u  field  in  case 
IV  after  three,  six,  nine  and  twelve  hours  of  integration. 

Fig.  19  is  also  representative  of  the  relative  performance  of  the  models  for  all  cases 
in  the  rotating  system.  Not  suprisingly,  the  results  of  these  comparisons  are  much  the 
same  as  in  the  non-rotating  cases.  SLSI  agrees  most  closely  with  the  high  resolution 
model.  FDEX  is  unable  to  capture  all  the  details  in  the  disturbance  fields  and  has  a 
tendency  to  place  the  extrema  slightly  upstream  from  the  HRSL  model.  The  SLEX 
forecasts  are  the  poorest  due  to  the  dominance  of  interpolation  errors. 
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Figure  7.  Experiment  1,  Case  I:  Total  u  Field  (nils)  after  A)  Two  B)  Four  C)  Six 
D)  Eight  Hours  of  Integration  for  HRSL  (solid  line),  SLS1  (dashed  line),  SLEX 
(broken  line)  and  FDEX  (dotted  line). 
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Figure  8.  Experiment  1,  Case  I:  </>'  Field  (nils)2  after  A)  Two  B)  Four  C)  Six  D) 
Eight  Hours  of  Integration  for  IIRSL  (solid  line),  SLSI  (dashed  line),  SLEX  (broken 
line)  and  FDEX  (dotted  line). 
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Figure  9.  Experiment  1,  Case  II:  Total  u  Field  (m/s)  after  A)  Two  B)  Four  C)  Six 
D)  Eight  Hours  of  Integration  for  HRSL  (solid  line),  SLSI  (dashed  line),  SLEX 
(broken  line)  and  FDEX  (dotted  line). 


25 


Figure  10.  Experiment  1,  Case  II:  </>'  Field  (mjs)1  after  A)  Two  B)  Four  C)  Six 
D)  Eight  Hours  of  Integration  for  HRSL  (solid  line),  SLSI  (dashed  line),  SLEX 
(broken  line)  and  FDEX  (dotted  line). 
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Figure  11.  Experiment  1,  Case  III:  Total  u  Field  ( m/s )  after  A)  Two  B)  Four  C) 
Six  D)  Eight  Hours  of  Integration  for  HRSL  (solid  line),  SLS1  (dashed  line),  SLEX 
(broken  line)  and  FDEX  (dotted  line). 
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Figure  12.  Experiment  1,  Case  III:  4>'  Field  (m/s)J  after  A)  Two  B)  Four  C)  Six 
D)  Eight  Hours  of  Integration  for  HRSL  (solid  line),  SLS1  (dashed  line),  SLEX 
(broken  line)  and  FDEX  (dotted  line). 
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Figure  13.  Experiment  1,  Case  IV:  Total  u  Field  (m/s)  after  A)  Two  B)  Four  C) 
Six  D)  Eight  Hours  of  Integration  for  HRSL  (solid  line),  SLSI  (dashed  line),  SLEX 
(broken  line)  and  FDEX  (dotted  line). 
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Figure  15.  Experiment  1,  Case  V:  Total  u  Field  (m/s)  after  A)  Two  B)  Four  C) 
Six  D)  Eight  Hours  of  Integration  for  HRSL  (solid  line),  SLSI  (dashed  line),  SLEX 
(broken  line)  and  FDEX  (dotted  line). 
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Figure  16.  Experiment  1,  Case  V:  4>'  Field  (m/s)1  after  A)  Two  B)  Four  C)  Six 
D)  Eight  Hours  of  Integration  for  HRSL  (solid  line),  SLS1  (dashed  line),  SLEX 
(broken  line)  and  FDEX  (dotted  line). 
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Figure  17.  experiment  2,  Case  I:  u  and  /  ( m/s )  after  A)  Four  B)  Eight  Hours  of 
Integration  with  1IR.SL. 
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Figure  18.  Comparison  of  Perturbations  in  Rotating  and  Non-rotating  System  (m/s)2 
for  Case  II  A)  u  and  B)  </>'  with  HRSL.  Solid  line  indicates  non-rotating  system;  the 
dashed  line  indicates  solution  for  the  rotating  system. 
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Figure  19.  Experiment  2,  Case  IV:  Total  u  Field  (m/s)  after  A)  Two  B)  Four  C) 
Six  D)  Eight  Hours  of  Integration  for  HRSL  (solid  line),  SLS1  (dashed  line),  SLEX 
(broken  line)  and  FDEX  (dotted  line). 
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V.  CONCLUSIONS 


This  investigation  is  a  direct  comparison  of  a  semi-Lagrangian,  semi-implicit  model 
to  the  semi-Lagrangian,  explicit  and  leapfrog  finite  difference  models  for  the  shallow 
water  system  with  bottom  topography.  Of  particular  interest  is  the  models'  ability  to 
simulate  hydrolaulic  jumps. 

In  the  first  experiment,  the  Coriolis  parameter  is  set  to  zero,  and  the  cases  which 
Petroliagis  (1988)  exai. lined  with  a  Galerkin  finite  element  model  are  reconstructed,  al¬ 
beit  on  a  larger  domain.  Initially,  the  How  was  uniform  and  zonal,  and  the  free  surface 
was  fiat.  According  to  the  theory  of  Houghton  and  Kasahara  (I96S),  the  first  two  cases 
are  subcritical.  A  stable  speed  maximum  should  appear  over  the  ridge.  The  last  case  is 
supercritical,  and  a  speed  minimum  should  develop  over  the  ridge.  The  remaining  cases 
are  unstable  and  should  produce  a  jump  in  the  vicinity  of  the  obstacle. 

The  SLS1  forecasts  are  superior  to  the  FDEX  and  SLLX  forecasts  for  all  cases. 
SLLX  forecasts  are  destroyed  by  excessive  smoothing  due  to  interpolations;  FDEX 
breaks  down  when  the  non-linear  interaction  become  too  large.  SESI  results,  however, 
are  consistent  with  hydraulic  jump  theory.  In  the  no  jump  cases,  it  forecasts  stable  speed 
and  pressure  maxima  and  minima  where  they  arc  predicted  by  theory.  S LSI  forecasts 
increasing  winds  and  dropping  pressure  for  the  jump  cases.  However,  the  details  of  the 
hydraulic  jump  which  forms  in  case  III  are  only  apparent  when  the  resolution  is  in¬ 
creased,  and  in  Case  IV,  the  domain  is  too  small  for  a  hydraulic  jump  to  form  before 
interference  from  the  transient  solutions  due  to  the  cyclic  boundary  conditions  occurs. 
Further  testing  should  be  done  on  a  larger  domain  (more  isolated  ridge)  with  the  higher 
resolution  semi-Lagrangian,  semi-implicit  model.  IIRSL. 

In  the  second  experiment,  the  same  cases  were  re-examined,  but  the  Coriolis  pa¬ 
rameter  for  45°N  was  used  to  examine  the  effect  of  rotation  on  jump  formation.  CRarlv, 
rotation  slightly  suppresses  the  perturbations  in  the  free  surface  and  delays  the  1  rma- 
tion  of  hydraulic  jumps,  but  further  study  is  necessary  for  a  more  quantitative  analysis. 
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